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for  an  oscillating  airfoil  and  for  the  transient  behavior  on  an^  g 

airfoil .where  a  large  separation  region  develops  or  when  laminar  flowlon 

breakdown  occurs.  The  algorithm  has  also  been  examined  for  three - - - - 

n/_ _ 

dimensional  marching  with  boundary  layer  and  Reduced  Navi er-S tokes  ty  Codes 

Av.iij.  and/or 

equations.  The  results  to  date  have  been  very  encouraging.  |Dist  Special 

i  I 


AzL 


AIR  FORCE  OFFICE  OF  SCIENTIFIC  RESEU^w  - 

NOTICE  of  transmittal  to  otic  {AFSC) 

This  technical  reoort  has  been  reviewed  an* 
approved  for  public  release  IAWAFR  190-12 
Distribution  is  unlimited.  2’ 

V  ' TTHF'.V  j.  ker-er 

’fctef.  Technical  Information  Division 

3 

j, 

v,  v. -V/.  v.  /•.  v.  a  v.  .  •’ .  r* 


1 

%\ 


as 


& 

■ft 


A.  STATEMENT  OF  WORK  (FROM  PROPOSAL)  2/1/85  -  10/31/86 

a)  A  composite  velocity/reduced  Navier-Stokes  (CV/RNS)  code  will  be 
developed  for  investigation  of  flow  past  lifting  airfoils  in  both 
subcritical  and  transonic  flow.  Full  viscous- inviscid  interactions  shall  be 
evaluated  for  laminar  and  turbulent  flow  conditions.  The  CV/Euler  solutions 
shall  be  critically  evaluated  for  accuracy,  prediction  of  entropy  and 
vorticity  patterns  and  variances  from  the  full  potential  solutions. 
Comparisons  will  be  made  with  the  CV/RNS  solutions  to  assess  the  effect  of 
viscous  interaction  on  the  inviscid  solutions. 

b)  A  time-consistent  relaxation  procedure  will  be  developed  for  modeling 
transient  behavior.  The  procedure  will  be  tested  on  airfoil  geometries  with 
transient  boundary  conditions  or  for  cases  where  transient  behavior  occurs 
due  to  large  Reynolds  number  effects. 

B.  STATUS  OF  RESEARCH  EFFORT  2/1/85  -  2/1/86 

As  given  in  the  Work  Statement  the  goals  and  status  to  date  can  be 
summarized  as  follows: 

(1)  Application  of  the  Composite  Velocity  (CV)  formulation  for  the 
Euler  equations  and  comparison  with  the  full  potential  solutions.  This  has 
essentially  reached  completion,  see  write-up,  B.l,  enclosed. 

(2)  Composite  Velocity/Reduced  Navier-Stokes  solutions.  This  has  been 
completed  for  steady  transonic  flow  over  a  NACA  0012  at  zero  incidence,  see 
reference  C.3  and  reference  E.2. 

(3)  Application  of  the  CV  formulation  to  lifting  airfoils.  The 
formulations  of  (1)  and  (2)  await  the  development  of  a  grid  generation 
technique  that  is  currently  under  investigation. 


(4)  Development  of  a  time-consistent  CV  relaxation  procedure  for  two- 
dimensional  transient  flow.  This  has  been  completed,  see  write-up  B.2  and 
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.  B.  3.  Preliminary  tests  have  also  been  considered  for  three-dimensional 

space  marching  (a  task  not  scheduled  until  the  option  period  11/1/86  - 
10/31/87),  see  write-up  B.2. 

(5)  A  new  iterative  'direct-solver'  formulation  to  be  applied  for  more 
efficient  computation  has  also  been  developed  and  is  discussed  in  B.4. 

Generally  speaking,  research  progress  has  been  quite  good  and  initial 
solutions  for  all  items  in  the  Work  Statement  should  be  completed  by 
10/31 /86. 


B. 1  COMPOSITE  VELOCITY  EULER/POTENTIAL  SOLUTIONS 


Solutions  to  the  Reduced  Navier-Stokes  equations  obtained  with  the 
composite  velocity  system  were  presented  at  the  22nd  Aerospace  Sciences 
Meeting,  January  1986,  see  reference  E.2.  It  was  noted  that  with  the 
formulation  of  that  time  entropy  was  generated  only  by  viscoscity  and  that 
entropy  was  not  generated  by  shock  waves.  The  following  describes  a 
procedure  for  obtaining-  Euler  solutions  for  the  outer  (Euler)  flow  with 
accurate  entropy  generation  due  to  shock  waves. 

The  composite  velocity  representation  may  be  written  as 


u  -  (1+<j>  ) 
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The  interpretation  of  the  composite  velocity  terms  for  inviscid  flows  varies 


slightly  from  that  for  viscous  flows.  The  $  term  still  represents  an 
irrotational  "pseudo"  potential  function.  The  U  term,  however,  is  no  longer 
associated  with  the  viscous  effects  but  rather  it  is  associated  with  the 


effects  of  rotational ity,  i.e.  vorticity  C  »  f(U,U  ). 


The  following  set  of  Euler  equations  are  obtained  for  2-D  flow 
Continuity 
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where  S  is  entropy  and  G  is  a  Bernoulli  like  pressure  term.  All  5“ 
derivatives  are  differenced  using  backward  differences  except  for  4^  which 
is  forward  differenced.  In  this  manner  4>  is  being  relaxed;  U  is  marched. 
This  differencing  is  consistent  with  Helmholtz's  vorticity  theorem. 

For  transonic  flows  the  Enquist-Osher  flux  biasing  scheme  has  been 
adapted  for  the  composite  velocity  system.  Thi3  scheme  consists  of  defining 
a  new  density 
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(pq)_  -  0.0  MSI,  (pq)_  -  pq  -  p  q  M  >  1. 

*  * 

Here  p  and  q  are  the  sonic  velocity  and  density.  This  procedure  produces 
very  sharp  shocks  and  guarantees  that  no  expansion  shocks  can  occur. 
Density  shifting  is  performed  only  on  terms  without  U  in  the  5_deri vatives 
in  the  continuity  equation. 

Equations  2-4  are  the  full  Euler  equations;  the  transonic  solutions 
calculated  using  these  equations  should  contain  the  rotational  effects  and 
the  correct  entropy  rise  at  the  shock  wave.  If  the  equations  are  solved  in 
their  present  form,  however,  no  entropy  is  generated  at  the  shock  wave  and 
the  i3entropic,  irrotational  full  potential  solution  is  obtained.  If,  on 


the  other  hand,  vorticity  is  introduced  into  the  flow  being  solved,  this 


additional  generation  of  vortlclty. 

The  reason  this  procedure  falls  to  capture  rotational  effects  lies  in 
the  form  of  the  5-momentum  Equation  3.  In  order  to  more  adequately  capture 
the  (Euler)  shock  wave  in  transonic  flows,  the  ^-momentum  equation  is 
rewritten  in  a  quasi  conservation  form 


p  +  £  [(ph2(U2+2U)uJ)£  ♦  (Ph1Uuev)n]  ♦  1  PuevUh1 
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For  transonic  flow,  values  of  U  will  be  generated  at  the  shock  wave.  Also, 
the  generation  of  U  values  will  bring  about  a  corresponding  entropy  rise  at 
the  shock  wave. 

The  modifications  to  the  5-momentum  equation  are  easy  to  impliment  in 

the  numerical  solution  procedure.  The  right  hand  side  of  Equation  3  is 

evaluated  at  the  previous  iteration  level.  The  (p  +  pu2)  term  may  be 

evaluated  using  a  two  point  backward  difference  and  a  second  order  central 

difference  is  used  for  the  (pu  v)  term.  All  the  other  terms  are  evaluated 

e  n 

at  the  (i,j)  location. 

While  this  form  of  the  equations  produces  the  required  entropy  rise  at 
the  shock  wave,  it  also  produces  spurious  entropy  in  the  regions  of  high 
gradients  at  the  leading  and  trailing  edges  of  the  airfoil.  Recalling  that 
the  nonconservation  form  produces  no  entropy,  but  only  convects  vorticity 
introduced  into  the  system,  this  form  of  the  equations  will  be  used 
everywhere  except  in  the  shock  region,  where  the  quasi  conservation  form  of 
the  equations  is  used.  The  shock  region  is  defined  as  the  region  from  the 


peak  Mach  number  to  the  point  where  the  change  in  pressure  is  less  than  one 


characteristic  that  no  entropy  Is  produced  except  in  the  shock  wave  region 
and  this  entropy  Is  convected  accurately  downstream. 

One  result  for  the  pressure  coefficient  along  a  NACA0012  airfoil  for  a 
freestream  Mach  number  of  0.85  is  presented  in  Figure  B1.1.  As  expected  the 
Euler  shock  is  weaker  and  lies  upstream  of  the  potential  shock.  The  Mach 
number  ahead  of  the  shock  has  a  value  of  1  .3395.  The  calculated  value  of 
the  Mach  number  after  the  shock  is  0.7683  and  the  entropy  generated  has  a 
value  of  S-0.0270.  These  values  compare  well  with  the  values  H^0.766H  and 
S-0.0283  obtained  from  the  shock  jump  conditions,  see  Figure  B1 .2.  There  is 
a  small  jump  in  the  value  of  entropy  at  the  trailing  edge  and  in  the  shock 
structure.  However,  the  former  effect  is  localized  to  the  trailing  edge 
region  and  arises  due  to  the  singularity  in  the  present  transformation. 
The  latter  overshoot  is  physical  and  typical  of'the  entropy  behavior  in  the 
shock  structure.  This  procedure  is  now  being  applied  for  three-dimensional 
and  other  two-dimensional  geometries  for  the  reduced  Navler-Stokes 


equations. 


B. 2  CONSISTENT  STRONGLY  IMPLICIT  ITERATIVE  PROCEDURES  FOR  TWO-DIMENSIONAL 
UNSTEADY  AND  THREE-DIMENSIONAL  SPACE-MARCHING  FLOW  CALCULATIONS 

P.K.  Khosla  and  S.G.  Rubin 
(To  be  submitted  for  publication) 

Introduction 

The  success  of  a  numerical  formulation  for  the  solution  of  two- 
dimensional  time  accurate  or  three-dimensional  spatial  accurate  flow,  e.g., 
RNS  marching  or  global  relaxation,  depends  largely  on  the  solution 
algorithm.  A  non-iterative,  unconditionally  stable,  consistent  procedure 
should  provide  maximum  efficiency.  For  large  time  increments  (At)  i.e., 
steady  state  calculations,  consistency  is  not  critical;  rather,  the 
technique  should  have  strong  convergence  properties.  For  moderate  At  and 
transient  flows,  or  moderate  Ax  for  spatial  marching,  consistency  of  the 
numerical  scheme  plays  an  important  role. 

The  most  generally  applied  implicit  and  consistent  formulation  is  the 

ADI  factorization  due  to  Douglas  and  Gunn^.  The  Briley/McDonald^  and 

Beam/ Warming^  schemes  are  based  on  this  ADI  factorization.  There  are, 

however,  a  number  of  problems  with  the  ADI  technique  that  have  been 

encountered  by  various  investigators.  These  are  (i)  an  "instability" 

associated  with  the  boundary  condition,  which  may  be  attributable  to  the 

choice  of  intermediate  boundary  conditions  and  can  only  be  controlled  with 

smaller  values  of  At,  (ii)  poor  rates  of  convergence  for  At  >  1,  i.e.  steady 

(4)  ( c)  2 

state  problems  (see  Nietubicz  and  Ghia  et  al  )  and  (iii)  0(At  ) 

accuracy  of  the  factorization,  which  may  not  be  desirable  with  certain  types 


of  discontinuous  initial  conditions.  For  the  ADI  method  and  also  Crank- 
Nicholson  scheme,  it  can  be  shown  that  odd  and  even  time  steps  decouple  for 
larger  At.  Thus  the  solution  is  oscillatory.  Additional  time  smoothing  is 
required  in  such  cases  to  make  the  technique  stable  and  convergent.  For 
At  >>  1,  inconsistent  and  0(At)  techniques,  e.g.  CSIP,  typically  possess 

(5) 

better  convergence  properties.  Ghia  et  al.  have  presented  comparisons  of 
the  CSIP  and  the  CADI  method  for  both  steady  and  unsteady  internal  flows. 
The  advantages  of  the  CSIP  for  steady  calculations  has  been  reaffirmed. 
Although  the  CSIP  has  also  been  used  by  these  investigators  for  unsteady 
flow  computations,  it  may  be  shown  that  the  consistency  and  hence  the  time 
accuracy  of  the  CSIP  will  deteriorate  as  the  spatial  mesh  width  h  is 

refined,  i.e.,  for  a  fixed  time  step,  >>  1  .  A  similar  limitation  exists 

h 

for  other  relaxation  procedures,  e.g.  LSOR,  SSOR.  This  deficiency  of 
relaxation  techniques  has  not  been  discussed  in  any  of  the  earlier  time 
accurate  references  cited  here.  Therefore,  there  is  currently  no  single 
implicit  solution  algorithm  that  has  both  the  desirable  properties  of 
consistency  (small  At),  rapid  convergence  to  the  steady  state  (At  >>  1)  and 
can  also  render  an  0(At)  scheme  consistent  so  as  to  be  applicable  when  first 
order  accuracy  is  desirable. 

In  the  present  investigation  a  number  of  simple  remedies  have  been 
investigated  to  render  any  of  the  well  known  relaxation  procedures 
consistent.  These  are  (i)  a  modified  predictor-corrector  SIP  or  C„^P 
procedure  (ii)  a  multi-grid  predictor-corrector  scheme  and  (iii)  a  new 

algorithm  based  on  the  Sherman-Morrison^  formula.  The  first  two  remedies 
require  considerable  programming  steps  and  results  in  two  step  procedures 
similar  to  the  ADI  technique.  The  third  formulation  is  very  simple  to 


implement,  requires  little  modification  of  existing  codes  and  results  in 
less  than  5%  additional  computational  effort.  Furthermore,  it  is  a  single 

2 

step  procedure  that  can  be  applied  to  to  achieve  either  O(At)  or  0(At  ) 
accuracy.  Intermediate  boundary  conditions  are  not  required  for  either  of 
the  procedures.  All  of  the  algorithms  have  been  tested  successfully  with 
the  SIP  on  a  model  problem.  These  algorithms  have  been  investigated  for 
application  to  flow  over  an  expanding  airfoil. 

MODIFIED  SIP  OR  CSIP 

The  new  procedures  are  described  here  for  the  model  problem: 

$  »  vV^<j>  (1) 

A  second  order  accurate  C rank-N ichol son  scheme  (6  =  1/2)  leads  to  the 
following  discrete  form  of  the  equations: 

— -  vA(e$n+1  *  ( 1 -9 ) qn) i j ,  (2) 

At 

where  the  various  diagonals  in  the  matrix  operator  A  are  typically  of 
2 

0(1 /h  );  h  is  the  mesh  spacing  and  n  the  temporal  index.  This  equation  can 
be  rewritten  as: 

(I  -  vAto  A)jn  »  (I  ♦  Atv(  1 -5 )  A)  J>n  ,  (3) 

n+i 

An  exact  solution  for  t  is  obtained  when  the  coefficient  matrix  ( I-vAt9A) 
can  be  inverted.  For  matrices  arising  from  general  two  dimensional  (Navier- 
Stokes,  Euler  or  RMS)  operators,  this  requires  the  use  of  Gaussian 
elimination  or  a  variant  thereof.  For  large  systems,  where  many  length 
scales  must  be  resolved,  such  a  solution  algorithm  is  extremely  inefficient 
and  prohibitively  expensive.  An  approximate  factorization  is  usually 
preferable.  This  factorization  should  be  efficient  and  ideally  generate  an 


approximate  solution  which  is  within  the  truncation  error  of  the  scheme. 


The  CSIP  is  such  a  procedure.  Application  of  the  SIP  for  the  inversion  of 


(I  -  vAt0A)  introduces  the  diagonal  elements  or  "corner  points" 


.  ,«  etc-  These  are  treated  explicitly  and  iterated  upon  in  order  to 

achieve  the  converged  solution.  A  closer  examination  of  the  coefficients  of 
these  terms  in  the  inversion  algorithm  reveals,  that  as  noted  previously 

these  diagonals  are  of  0(^~) .  Thus  an  accurate  and  consistent  prediction 

h 


of  <}>?+1  requires  that  <<  1. 

h 


This  restricts  the  choice  of  At  to 


unacceptable  small  values  for  h  <<  1.  In  the  following  section,  a  number  of 
simple  remedies  have  been  investigated  for  improving  the  consistency  of  the 
SIP  procedure.  These  remedies  are  quite  general  and  can  be  used  with  other 
relaxation  procedures  too. 


Predictor-Corrector 

In  order  to  apply  a  single  solution  algorithm  for  both  steady  and 
unsteady  or  3“D  marching  problems,  i.e.,  a  single  formulation  retaining  the 
most  desirable  features  for  all  values  of  At  (or  Ax  marching),  the  present 
authors  have  investigated  a  number  of  two-step  and  one-step  techniques. 
Intermediate  boundary  conditions  are  not  required. 

(i)  The  matrix  A  resulting  from  the  spatial  or  "non-marching" 
discretization  is  rewritten  as: 

A <J>  »  vU"1  +  1  +  ( 1  — e ) )  =  M<j>n  1  +  ( 1  —  e ) N4>n .  (4) 

yy  z  z  z  z 

Equation  (2)  can  then  be  written  as: 

(I-vAt0MHn+1  -  (I  +  vAt(  1  -9 )  A)  <4>n  +  9(  I-e  )N<f>n  (5) 
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In  the  predictor  step,  the  coefficient  matrix  (I  -  vAt9M)  with  e 


e  and 
o 


9  -  1  is  Inverted  by  the  SIP. 


The  error  arising  from  the  diagonals  $ 


n 

1-1 ,J*1 


and  01;  can  be  made  independent  of  the  mesh  size  by  an  appropriate 

i  *  *  ,  J  “  ' 


.  n*  1 


choice  of  the  parameter  e.  This  provides  a  reasonable  predictor  for  ♦ 

The  error  due  to  the  inconsistency  is  within  the  truncation  error  of  the 
numerical  approximation  of  the  scheme.  A  second-order  accurate  corrector 


step  repeats  the  procedure  with  eq  -  1  and  0  -  The  choice  of  the 


appropriate  value  of  eq  is  crucial  to  the  success  of  the  method.  The 
following  expression  was  selected  in  the  test  problem: 


At/h2 

1  ♦  At/h2 
Other  expressions  or  values  of  e 

o 


(6) 

can  also  be  chosen  and  still  provide  the 


2 

required  consistency.  It  should  be  noted  that  for  (At/h  )  <<  1,  the  order 
of  inconsistency  is  bounded  and  independent  of  the  mesh  spacing  h.  For 
2 

At/h  >>  1 ,  e  -  1  and  the  usual  SIP  algorithm  is  recovered.  A  model  problem 
with  Dirichlet  boundary  conditions  was  chosen  to  test  the  predictor- 

corrector  SIP  technique.  Calculations  for  v  •  —  were  considered  for 

(17x17)  and  (51x51)  grids  and  At  -  1.  For  a  single  time  step,  the  maximum 
error  in  the  numerical  solution  is  less  than  2\  and  71,  respectively. 
Additional  corrector  iterations  with  e  »  1  will  further  reduce  this  error. 
For  comparison,  the  standard  SIP  on  a  (51x51)  grid  incurs  almost  a  1 00% 
error .  The  method  requires  the  splitting  of  various  derivatives  in  a  manner 


similar  to  the  ADI  method.  The  choice  of  e  is  quite  arbitrary  and  may  even 
depend  upon  the  nature  of  the  differential  equation  and  the  flow  parameters. 

SIP-Multlgrld 

A  complementary  approach  to  improve  the  consistency  of  the  CSIP 
involves  the  application  of  a  multigrid  procedure.  Since  consistency  error 
is  amplified  for  fixed  At  and  fine  meshes,  a  course  mesh  predictor  minimizes 
this  error.  The  predictor  step  is  then  repeated  on  a  succession  of  finer 
grids,  the  fine  grid  corrector  step  will  not  require  additional  iterations. 

2  2  2 

Solutions  obtained  in  this  fashion  are  also  0(At  ,Ay  ,Az  );  however,  local 
convergence  is  enhanced  by  the  multigrid  smoothing  process.  Results  similar 
to  the  one  obtained  in  the  previous  method  are  also  obtained  by  this 
technique.  However,  the  programming  complexity  increases.  This  technique 
can  always  be  employed  with  the  procedures  described  in  this  paper.  A  new 
iterative  procedure  that  is  based  on  this  idea  and  the  application  of  sparse 
matrix  direct  solver  is  being  currently  investigated  (reference  C4)  for  the 
solution  of  the  reduced  Navier-Stokes  equations.  The  preliminary  results 
are  quite  encouraging. 

C .  Stone*  s-SIP 

The  consistency  of  the  SIP  procedure  can  be  improved  considerably  by 
using  a  second-order  factorization  similar  to  the  one  originally  proposed  by 
Stone.  In  the  present  section  a  number  of  techniques  similar  in  character 
to  Stone's  procedure  will  be  discussed.  Although  most  of  these  techniques 
are  unstable  for  steady  state  calculations  (some  more  than  the  others),  they 
provide  simple  extensions  of  SIP  for  time  consistent  computations,  fairly 


large  At  -  1  and  highly  stretched  grids.  Some  additional  inexpensive 


remedies  are  suggested  which  not  only  improve  the  stability  for  larger  At, 
but  also  improve  the  consistency. 


Stone  had  proposed,  the  following  second-order  factorization.  This 
cancels  the  corner  values  of 

<tn+^  -  <fcn  —  ct<  <bn  +  <bn  -  <hn  ) 

Vl,j+1  i-1 ,  j+1  avVl,j  *i,J+1  ’ij' 

♦  aUn+1  *  d,n+1  -*n+1) 

vVl,j  ^i,  j  +  1  *i,j' 


cn+1 

Pi+1,j-1 


,j-1  "  o(*i+1,J  +  *i.J-1 


,.n+1  A  ,n+1  .n+1 . 

*  *  ♦i.j.i  -  *i,j> 


Clearly  the  terms  evaluated  at  the  previous  time  level  have  a  spatial 
truncation  error  of  O(AxAy)  when  a  is  chosen  to  be  unity.  This  happens  to 
be  the  simplest  method  for  making  the  method  consistent.  Originally, 
Stone's  technique  was  abandoned  because  of  the  complexity  of  the  factorizing 
procedure  for  more  than  one  unknown..  However,  the  following  simple 
implementation  has  been  found  to  be  equivalent. 


Given  the  algebraic  system 

Mi.j-l  +  Vi-I.j  +  B1*i,j  +  C1  ^i,  j+1  +  E1*i+1,J  (7) 

the  solution  algorithm  can  be  described  as: 

*i.J  *  +  Ei. j^i. j+1  +  Fi,jVl,j  * 

The  elimination  of  the  lower  triangular  terms  .  .  and  .  along  with 

1  '  »  J  *  *  J“  * 

the  cancellation  of  the  corner  points  <t>.  .  and  4> .  .  can  be  carried 

•  9  j*“ i  i*  i  f  J  + 1 

out  in  a  single  step  as: 
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% 

& 


1 

♦ 

va 

I 


♦iij-1  ‘  (1  '  aFi,j-ir1  [GMi,j-1  +  '  aFi,j-1Hi,j 

+  aFi,j-1  Vl,j  +  ^i+1 , j-1  '  fl,(*i+1fJ^i,J-r*lfJ))J 

(8) 

and  a  similar  expression  for  4»i_1  j.  Substituting  these  in  the  governing 

equations,  we  get,  the  following  recurrence  relation  for  F.  (say), 

1  9  J 

Fl’J  B  •  4 

’  '  ’  -  '  ’  -  aEl-,,J 

It  may  be  noticed  that  the  values  of  F.  .  do  not  depend  upon  the  grid 

i » J 

spacing  and  F  /(1-F.  )  etc.  are  0(1).  The  group  of  terms  being 

l ,  J  1  • J“ i 

computed  explicitly  in  equation  (8)  are  thus  of  O(AxAy).  As  such,  the  error 
does  not  grow  in  a  space  maraching  or  time  dependent  problems  due  to  the 
initial  guess  of  the  solution.  A  general  block  version  of  this  algorithm 
for  any  number  of  equations  and  unknowns  has  been  coded  and  as  a  test  is 
being  applied  to  the  three  dimensional  boundary  region  equations  for 
supersonic  flow  past  a  cone.  A  second  technique  allows  for  the  inclusion  of 
an  additional  diagonal  (corresponding  to  #1+1  j  +  1)  in  the  sparse  LU 

factorization  of  the  preconditioning  matrix  M.  In  this  case,  the  source 
terms  computed  at  the  nth  time  level  can  further  be  cancelled  by  addition 
and  subtraction  of  a  grouping  of  the  form  [$.  .  .  .  -  a(<t>.  .  ..." 

<j>.  ,)].  The  storage  required  for  the  factorization  step  in  this  case 
*  >  J 

increases  from  3N  to  4N.  Other  groupings  can  be  devised  to  achieve  similar 
results.  For  a  -  0,  all  the  methods  reduce  to  the  simple  SIP  which  has  been 


V.  'V 


used  by  the  present  authors  for  coupled  system  and  for  a  large  number  of 


steady  state  computations. 

The  cancellation  of  the  corner  points  ♦,  .  .  .  and  ,  .  .  can  also  be 

I  »  J4,  1  9  J""  1 

performed  by  changing  the  truncation  error  of  the  governing  equation.  For 
example,  the  conservation  form  of  the  following  equation 

0 


+  A  +  B 
3t  x  y 


where  A  and  B  are  functions  of  #,  <p  and  0  ,  etc. 

x  y 


Then, 

(Ax)j.,  .  (By),,,  -  (ixljj  -  (Byjjj 
(Ax)j„  .  (By)1H  -  (AxJj  j  -  < By > j ( j 

can  be  used  to  cancel  the  corner  values  of  j  +  1  and  j_1 .  This 

procedure  worked  quite  well  for  steady  subsonic  potential  flows.  Additional 
iterations  with  a  -  1  did  not  diverge  when  used  to  compute  flow  past  a 
biconvex  airfoil.  However,  this  may  be  fortuitous.  Unfortunately,  this 
technique  is  quite  complicated  and  becomes  prohibitively  cumbersome  for  more 
than  one  unknown. 

In  order  to  eliminate  the  sensitivity  of  Stone’s  procedure  to  the  value 
of  a,  a  rank  one  improvement  of  the  iterative  procedure  was  found  to  be 
quite  useful.  This  is  achieved  by  the  application  of  a  Sherman-Morrison 
formula.  It  must  be  emphasized  that  all  iterations,  required  at  a  given 
time  step  for  nonlinear  convergence,  can  be  performed  with  smaller  values  of 
a.  A  value  of  a  >  0.9  has  been  utilized  for  this  purpose. 
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Sherman-Morrison  Formula 


This  formula  inverts  a  matrix  of  the  following  form: 

A  -  B  +  (JVT 

The  matrix  B  is  such  that  B  1  can  easily  be  computed  and  U  and  V  are  two 
vectors.  Then 


-1  T 
B  UV  B 

T 

1  ♦  V  B 


-1 


U 


(9) 


> 


Wilfv  '  used  this  formula  to  invert  any  non-singular  square  matrix.  The 

computational  cost  is  O(N^)  which  is  comparable  to  Gaussian  elimination. 

2 

Memory  requirements  are  0(N  )  for  N  large.  As  a  direct  solver  the  Sherman- 
Morrison  Formula  is  not  competitive  with  iterative  techniques.  However,  it 
can  be  usefully  employed  for  improving  the  consistency  of  iterative 
techniques.  Most  Iterative  methods  are  based  on  some  type  of  splitting  of 
the  coefficient  matrix,  e.g., 

A$  ■  b  .  (10) 

can  be  written  as 

(M+NH  -  b  (11) 


I 


F 


-1 


For  most  of  the  matrices  arising  in  flow  problems,  M  can  be  computed  with 
reasonable  computational  effort;  however,  the  error  matrix  N  cannot  be 
decomposed  in  the  required  form,  i.e., 

,T 


N  »  UV  ‘  (12) 

Therefore  equation  (9)  cannot  be  usefully  exploited.  However,  the  linear 
equation  (11)  can  be  written  as: 

(M  ♦  N  <j><J>T)<i>  -  b  ,  (13) 
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vlvlvl  Ivlvlv.'vi:  . 


where  0  -  <p/<<j>  ,  o>  is  the  unit  vector.  Equation  (13)  is  non-linear  and  can 
be  solved  iteratively  for  $  as: 

(M  ♦  U$T  )4>n*1  -  b  (I1*) 

where 

A 

U  -  (NO) 

The  solution  of  (I1*)  can  be  written  as: 

(M_1 U ) <$TM~1  b> 

4.  -  M_1b - -  (15) 

1  +  <$TM  U> 

This  solution  further  improves  the  symmetry  property,  as  well  as,  the  time- 
consistency  of  the  solution.  From  equation  (15),  it  can  be  seen  that  the 
additional  effort  required  in  the  present  case,  as  compared  to  the 

relaxation  procedure  based  on  the  preconditioning  by  M  1 ,  is  associated  with 

the  evaluation  of  two  scalar  products.  M  ^b  and  M  ^  can  be  computed  in  the 
same  loop  with  very  small  increase  in  the  computation.  The  overall  increase 
in  time  is  between  2%  to  5% .  The  maximum  error  is  less  than  2%  for  a 
(17x17)  grid  with  At  -  1  .  The  problem  under  consideration  has  smooth 
initial  conditions,  with  the  exact  solution  having  a  0fc  of  0(1).  For  a 

problem  with  sudden  heating  of  the  boundary,  the  initial  conditions  are 
discontinuous.  In  such  a  case,  time  consistency  of  the  solution  requires 
the  use  of  a  smaller  At,  at  least  in  the  initial  stages  of  the  calculation. 


The  error  matrices  for  SIP  with  and  without  the  Sherman-Morrison  update  are 


No  and  (No)  respectively.  Since 
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the  influence  of  the  error  matrix  is  reduced  with  the  application  of  the 


Sherman- Morrison  formula. 


Error  Analysis 

The  first  order  implicit  formulation  of  the  time  dependent  equation  can 
be  written  as: 


[I  +  At(M+N)J*n+1  -  ♦h  +  at  b 

where  the  index  corresponds  to  a  time  step  and  (M+N)  is  the  coefficient 
matrix  arising  out  of  the  spatial  terms.  The  time  consistency  of  the 
solution  procedure  with  Sherman-Morrison  update,  can  be  investigated  by 
rewriting  the  lefthand  side  of  the  above  equation  as: 


a  .  it  h  *  it  n  ;n  *  At  (»  -  n  »n 


n+1 


where  <f  represents  the  normalized  unit  vector  defined  earlier. 

Using  Taylor  series  expansion,  it  can  be  shown  that  the  error  term  is 
given  as 


N  *n 


<4>  ,  , 

«  f2  <,n> ;  v~> *  v, 

n*1  n+1 


-  (N  iTn+1  +  N  $n+1  +  0(At2) 

Clearly,  the  error  term  is  first  order  in  time  provided  is  0(1)  and  N  is 


such  that  the  spatial  error  does  not  dominate  the  error  term  as  the  grid  is 
refined.  In  such  a  case,  the  limits  At  +  0  will  render  the  system 
consistent.  The  strongly  implicit  procedure  with  a  »  0,  which  has  been 
extensively  used  in  previous  calculations,  does  not  lead  to  a  consistent 


splitting  of  the  coefficient  matrix.  The  error  matrix  N,  amplifies  the 

2 

error  as  h  +  0,  leading  to  a  limitation  that  (At/h  )  <<  1  for  consistency. 
However,  for  a  -  1,  the  spatial  error  is  not  amplified  as  h  +  0.  Thus 
consistency  is  achieved  if  At  and  h  *  0  independently.  The  procedure  is 
also  applicable  to  second  order  time  accurate  technique. 


Applications 


The  technique  has  been  tested  on  a  variety  of  simple  model  problems. 
Both  linear  diffusion  and  nonlinear  Burger  equations  have  been  investigated 
to  test  the  validity  of  the  procedure.  In  addition,  the  unsteady  flow  past 
a  biconvex  airfoil  was  also  considered.  The  results  for  simple  diffusion 
are  described  below. 

The  heat  conduction  equation  is  given  as: 


♦t  *  h 

The  exact  solution  to  this  problem  is 

4>(t,x,y)  »  e  sinBx  sinBy  +  e  e  sinBy 


Exact  initial  and  boundary  conditions  have  been  imposed  on  the  numerical 
solution.  This  problem  was  chosen,  because  the  boundary  conditions  and  the 
solution  are  time  dependent.  The  maximum  residue  on  17x17  and  51x51  uniform 
grids  are  depicted  in  Figures  B2.1  and  B2.2  for  both  the  consistent  SIP  and 
the  inconsistent  SIP  procedures.  No  additional  iterations  have  been 
performed  for  either  case.  Both  techniques  include  the  Sherman-Morrison 
procedure.  This  improves  the  accuracy  of  the  standard  or  inconsistent  SIP 
method  as  well  as  the  new  procedure.  On  coarser  grids,  the  two  results  are 


similar;  however,  on  the  fine  grid  there  is  significant  difference.  This 
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reflects  the  severe  inconsistency  of  the  standard  SIP  procedure.  The  new 
formulation,  however,  retains  the  error  of  the  order  of  truncation  error. 
This  is  even  true  for  finer  grids  and  large  values  of  At. 

As  another  example  of  the  applicability  of  the  new  algorithm,  the 
composite  velocity  solution  past  a  biconvex  heaving  airfoil  at  *  0.65  has 

also  been  computed.  Computations  have  been  carried  out  for  number  of  time 
steps  on  a  75*33  non-uniform  grid.  In  these  calculations,  the  non-linearity 
has  been  treated  via  picard  iteration.  At  least  one  additional  iteration  is 
performed  after  time  stepping  and  before  the  density  is  updated.  In  all  the 
calculations,  only  three  density  updates  have  been  performed.  For 
comparison  purposes,  solutions  using  a  direct  solver  have  also  been 
obtained.  Three  iterations  on  density  are  also  performed  for  these  exact 
computations.  The  results  of  these  calculations  are  depicted  in  Figs.  B2. 
Figures  (B2.3)  shows  the  coefficient  of  pressure  at  different  normalized 
time  levels.  These  calculations  have  been  performed  with  At  -  0.1  and  0.5. 
The  results  are  very  sim’lar  at  similar  time  levels.  For  comparison  the 
pressure  time  history  at  the  mid  chord  has  been  depicted  in  Fig.  B2.4  along 
with  a  direct  solver  solution,  see  B.4.  The  two  solutions  are  in  excellent 
agreement. 

As  another  test  problem  the  laminar  boundary  region  equations  for  flow 

over  a  10°  cone  at  Ma  ■  2  has  also  been  considered.  The  boundary  region 

equations  contain  all  cross  flow  diffusion  terms  and  are  solved,  using  the 
coupled  algorithm.  Solution  for  this  case  is  well  known.  The  present 
consistent  CSIP  requires  only  two  iteration  as  compared  to  5  to  6  iteration 
with  the  old  CSIP  to  achieve  the  same  order  of  accuracy  for  each  marching 
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B.3  A  COUPLED  STRONGLY  IMPLICIT  PROCEDURE 
FOR  REDUCED  NAVIER-STOKES  EQUATIONS 

Some  solutions  for  the  unsteady,  compressible  Reduced  Navier-Stokes 
equations  obtained  using  the  Sherman  and  Morrison  technique  were  presented 
at  the  AIAA  24TH  AEROSPACE  SCIENCES  MEETING  ,  held  in  Reno,  Nevada.  This 
algorithm  is  suitable  only  for  a  H-type  grid.  Since,  there  are  some 
inherent  advantages  in  a  C-type  grid  for  many  flow  problems,  an  algorithm 
that  can  be  used  for  both  the  H-type  and  the  C-type  grid  would  be  more 
desirable.  Such  an  algorithm  has  been  developed.  This  algorithm  is  based 
on  the  coupled  strongly  implicit  (CSIP)  procedure  of  Stones. 
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The  quasi- linearized  form  of  the  governing  finite  difference  equations 
together  with  the  boundary  conditions  can  be  written  as 
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,  I.  .  are  known  3X3  coefficient  matrices  , 
1  f  J 


J.  is  a  known  3X1  matrix  and  '  n'  refers  to  the  present  time  level. 

1  t  J 

The  CSIP  is  developed  from  the  following  approximate  LU  decomposition. 
Vn  -  p  +  q  Vn  +  R  Vn 

i.j  'i.j  ui,j  vi,j-l  Ki,j  Vl.j  {2) 
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[vV  the  recursion  relation 


i ,  j  ~  -  '  '  i,J+1  *  Qi,J+1  ’  Ri,j+1  * 


Function  (  P. 


p . .  Qi*t,k  •  Ri.i.k  k-J  or  J*’  > 
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Similar  relations  are  obtained  for  Q.  .  and  R. 

*  »  J  +  9  J 


We  start  the  solution  procedure  by  solving  for  Pu  ,  ,  Q,.  .  and  R,,  .  . 

“ » J  M.J  M,j 


The  outflow  boundary  condition  (p  ).  ..  -  0  enables  us  to  solve  for  P„  . 

X  1*M  M,J 


Q„  .  and  R,,  .  from  the  governing  equations.  The  boundary  conditions  for  j-N 
M » J  M  v  j 


(p.  -1,u.-u  )  together  with  the  continuity  equation  for  j-N-1/2 
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P,  „  ,  Q.  „  and  R,  „  2  S  i  S  H  are  all  known  the  recursion  relation  can  be 
1  »  N  i  t  N  1 1  N 


used  to  obtain  P  ,  Q  .  and  R.  .  for  i-M-1 ,M-2, . . . . ,2  in  that  order.  The 
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continuity  equation  at  j-3/2  and  the  boundary  conditions  for  j-1 


(  p.  ,  -  1  ,  u.  ,  -  u  )  are  then  used  to  solve  for  v"  ,  .  From  p?  ,  , 
1,1  i,1<»  i,i  1,j 
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R.  .,2SjSN,2SiSM  one  can  solve  for  V  for  all  values  of  i  and 
1 » J  1 » J 


j  using  equation  (2). 

The  unsteady  breakdown  of  the  laminar  calculations  for  the  flow  past  a 
sine-wave  geometry  was  studied  using  this  algorithm  to  explore  the 
possibilities  of  obtaining  unsteady,  unsymmetric  solutions  for  large 
Reynolds  numbers.  So  far  such  a  solution  has  not  been  obtained.  Instead  it 
was  found  that  even  without  any  imposed  conditions  of  symmetry  in  the  wake 
the  solution  for  R-400,000  neither  converges  nor  exhibits  any  tendency 


towards  unsymmetric  shedding.  The  behaviour  in  time  of  the  wall  shear 
stress  for  R-400,000  is  shown  in  Fig.B3.l.  The  size  of  the  separated  region 


initially  increases  with  time  and  for  large  values  of  t  the  separated  bubble 
breaks  up  into  two  and  also  the  negative  peak  in  the  wall  shear  distribution 
increases  unboundedly  with  time.  Further  analysis  is  required  to  understand 
this  phenomenon  better. 

B.4  BLOCK  ITERATIVE  PROCEDURE  FOR  THE  SOLUTION  OF  FULL 
AND  REDUCED  NAVIER-STOKES  EQUATIONS 

Most  relaxation  methods  slow  down  when  the  number  of  grid  points 
increase.  Usually,  such  procedures  are  acceleralted  by  using  either 
conjugate  gradient  or  multi-grid  techniques.  For  algebraic  equations 
arising  from  large  Reynolds  number  flows  on  stretched  grids,  both  methods 
require  special  considerations  and  become  quite  problem  dependent.  There  is 
no  single  procedure  which  the  CFD  community  can  reliably  utilize,  without 
many  changes  in  the  codes,  that  will  work  for  a  large  class  of  problems. 
Usually  the  separated  flow  regions  require  special  attention  in  applying 
multi-grid  or  other  acceleration  techniques.  In  view  of  these  problems,  a 
new  technique  is  being  investigated.  This  is  based  on  a  direct  solver  and 
is  iterative  in  character.  The  relaxation  process  is  carried  out  by 
simultaneously  solving  the  equations  on  large  blocks  of  the  grid  and 
Iterating  between  these  blocks.  Typically,  direct  solvers  require  very 
large  amounts  of  memory  and  become  slow  as  the  numbe  of  grids  increases. 
However,  they  can  be  used  on  smaller  blocks  of  the  grid  in  a  fairly 
efficient  fashion.  In  this  manner  longer  wave  lengths  of  the  solution  error 
can  be  treated  more  effectively;  the  length  is  determined  by  block  size. 

The  block  iteration  procedure  consists  of  directly  solving  for  values 
cn  subgrids  of  the  main  grid.  It  has  been  found  experimentally  that 


selection  of  the  subgrids  and  the  order  in  which  the  iteration  takes  place 
has  a  significant  effect  on  convergence  rate.  The  method  has  been  tested  on 
the  Laplace  equation  in  the  unit  square.  The  fastest  grid  sequence  found  so 
far  is  described  in  Fig.  B4.1,  where  the  subgrid  numbers  show  the  sequence 
in  which  they  are  iterated.  The  first  two  subgrids  split  the  main  grid  in 
half  and  a  third  grid  that  overlaps  the  other  two  is  solved  last.  This 
subgrid  sequence  converges  in  6  iterations  (the  error  criterion  was  the 

difference  of  the  computed  and  exact  solutions  and  was  less  than  10  )  for  a 
21  by  21  grid,  and  7  iterations  for  a  101  by  101  grid.  Without  the  third 
grid  the  iteration  count  for  the  21  by  21  case  was  25. 

The  main  grid  was  further  divided  into  smaller  grids  as  shown  in  Fig. 
B4.2.  The  first  four  grids  were  solved  in  an  upward  sweep  followed  by  a 
downward  sweep  for  the  three  overlapping  grids.  This  was  done  to  ensure  a 
uniform  propagation  of  information  over  the  solution  domain.  As  expected 


) 


this  grid  strategy  needed  more  iterations  to  converge.  For  a  21  by  21  grid, 
15  iterations  were  required,  and  25  iterations  were  required  for  a  1 01  by 
101  grid.  The  convergence  histories  for  these  grid  strategy  in  a  model 
problem  are  shown  in  Figures  B4.3  and  B4.4  respectively.  Apparently  the 
method's  sensitivity  to  main  grid  size  increases  with  the  number  of 
subgrids. 


It  would  seem  initially  that  this  strategy  would  be  inferior  to 
directly  solving  for  the  entire  grid  (given  enough  computer  memory).  For 
nonlinear  problems  however,  iteration  is  required  anyway  to  converge  the 


nonl inearities.  The  amount  of  time  required  by  sparse  matrix  solvers 

1  5 

increases  superl inearily  with  the  number  of  unknowns  (0(n  ’  )  is  optimum). 


Therefore  one  iteration  of  this  subgrid  strategy  can  be  faster  than  an 
iteration  of  a  fully  direct  solution.  If  the  number  of  iterations  using  the 


subgrid  strategy  is  about  the  same  as  the  number  required  to  converge  the 
nonlinearities,  the  subgrid  strategy  can  be  faster.  The  procedure  has  been 
applied  to  the  solution  of  flow  in  a  driven  cavity  and  the  steady  laminar 
flow  past  a  NACA  0012  airfoil  at  Re  -  2000  and  =  .72.  Similar  rates  of 

convergence  have  been  obtained  for  these  problems. 

C.  TECHNICAL  PUBLICATIONS  2/1/85  -  2/1/86 

(1)  Ramakrishnan,  S.V.  and  Rubin,  S.G.,  "Numerical  Solution  of  Unsteady 
Compressible  Reduced  Navier-Stokes  Equations"  (to  be  submitted  to  AIAA 
Journal  or  Computers  &  Fluids). 

(2)  Khosla,  P.K.  and  Rubin,  S.G.,  "Consistent  Time  and  Space  Marching 
Coupled  Strongly  Implicit  Algorithms  "  (to  be  submitted  to  Computers 
and  Fluids  or  J.  Computational  Physics). 

(3)  Rubin,  S.G.  and  Gordnier,  R. ,  "Composite  Velocity,  Potential  Euler  and 
RNS  Solutions"  (in  preparation). 


(4)  Khosla,  P.K.  and  Bender,  E.,  "A  'Direct  Solver'  Iterative  Formulation 
and  Application  to  Viscous  Interacting  Flows"  (in  preparation). 


.  PROFESSIONAL  PERSONNEL  2/1/85  -  2/1/86 
S.G.  Rubin,  Professor  of  Aerospace  Engineering 
P.K.  Khosla,  Professor  of  Aerospace  Engineering 
S.V.  Ramakrishnan,  Research  Assistant  (PhD  student) 

R.  Gordnier,  Research  Assistant  (PhD  student) 

M.S.  "Transonic  Viscous  and  Euler  Solution  Using  a  Composite  Velocity 
Procedure",  March  1985. 

E.  Bender,  Research  Assistant  (PhD  student) 

E.  OTHER  ACTIVITIES  AND  INTERACTIONS  2/1/85  -  2/1/86 

( i)  Papers,  Presentations,  Seminars,  Short  Courses 

(1)  Ramakrishnan,  S.V.  and  Rubin,  S.G.,  "Numerical  Solution  of  Unsteady 
Compressible  Reduced  Navier-Stokes  Equations",  Paper  No.  AIAA  86-0205, 
AIAA  24th  Aerospace  Sciences  Meeting,  Reno,  NE,  January  1986. 

(2)  Gordnier,  R.,  "Transonic  Viscous  and  Inviscid  Solutions  Using  a 
Composite-Velocity  Procedure",  Paper  No.  AIAA  86-0074,  AIAA  24th 
Aerospace  Sciences  Meeting,  Reno,  NE,  January  1986. 

(3)  Reddy,  D.R.,  Delaney,  R.  and  Rubin,  S.G.,  "Reduced  Navier-Stokes 
Relaxation  Procedure  for  Three-Dimensional  Internal  Flows  with 
Interaction",  SAE  Aerospace  Technology  Conference,  Long  Beach,  CA, 
October  1985. 

(4)  Khosla,  P.K.  and  Rubin,  S.G.,  "Consistent  Strongly  Implicit  Iterative 
Procedures",  (submitted  to  10th  International  Conference  on  Numerical 
Methods  in  Fluid  Dynamics,  Beijing,  July  1986). 

(5)  Rubin,  S.G.,  "Global  Relaxation  Procedures  for  a  Reduced  Form  of  the 
Navier-Stokes  Equations" 


a.  Pennsylvania  State  University,  June  1985 


b.  Institute  for  Computer  Applications  in  Science  &  Engineering, 
(ICASE),  NASA  Langley  Research  Center,  Hampton,  VA,  May  1985 


c.  Allison  Gas  Turbine  Division,  General  Motors  Corporation, 
Indianapolis,  IN,  February  1985 

d.  University  of  Tennessee,  Knoxville,  TN,  July  1985.  Short  Course  on 
Finite  Element  Analysis  in  Fluid  Mechanics  and  Heat  Transfer. 

(6)  Khosla,  P.K.,  "Subsonic  and  Transonic  Reduced  Navier-Stokes  Techniques 

and  Applications",  International  Conference  on  Fluid  Dynamics,  Tokyo, 

Japan,  September,  1985. 

( ii)  Consulting  and  Advisory  Functions  2/1/35  ~  2/1/86 

S.G.  Rubin:  Principal  Investigator 

(a)  Joint  NASA/DOD  Panel  on  Hypersonic  Flow  Research  and  Graduate 
Training,  NASA  Headquarters,  May  1985. 

(b)  NASA  Lewis  Computational  Mechanics  Advisory  Committee.  To  advise 
the  Lewis  Research  Center  on  Computational  Mechanics  research  and 
development  of  a  Computational  Mechanics  Institute  under  the 
auspices  of  the  University  Space  Research  Association.  Meetings 
were  held  during  i 98U  and  1985. 

(c)  NASA  Aerospace  and  Research  Technology  Subcommittee  (ARTS)  of  the 
Office  of  Aeronautics  and  Space  Technology  (OAST).  Appointed 
December,  1985. 

(d)  Discussions  with  Dr.  J.  Shang  of  Wr ight-Patterson  Air  Force  Base 
for  interaction  on  CFD  and  for  sabbatical/consulting  arrangement 
starting  later  in  1986.  Agreement  in  principle  has  been  reached. 

(e)  Consultant  to  AVCO  Corporation,  Everett,  Mass  on  CFD  related 
problems,  November  1 98^  -  March  1985. 


(f)  Consultant  to  Allison  Gas  Turbine  Division,  General  Motors 
Corporation  on  use  of  Reduced  Navier-Stokes  Methodology  for 
problems  in  internal  flows,  October  1 98-4  -  Present. 

(g)  Advisor  to  Aerospace  Corporation  on  problems  associated  with 
viscous  interacting  flows  and  CFD,  September  1979  -  Present. 

F.  OTHER  INFORMATION 

The  Reduced  Navier  Stokes  (RNS),  Composite  Velocity  (CV)  and  Coupled 
Strongly  Implicit  Procedure  (CSIP)  have  been  and  continue  to  be  applied  for 
other  agencies  and  by  several  ether  investigators.  Recent  publications 
based  on  these  ideas  include:  R.H.  Pletcher  (Iowa  State  University)  for 
J.  Heat  Transfer  (to  appear  1986),  M.  Israeli  and  M.  Rosenfeld  at  the  AIAA 
6th  CFD  meeting  in  Cincinnati,  June  1985,  M.  Barnett  and  R.T.  Davis  in 
Computers  and  Fluids  (to  appear  1986),  H.  Raven  and  M.  Hoekstra  at  several 
hydrodynamics  symposia  (most  recently  in  Washington,  D.C.,  June  1985)  and  by 
B.  Laksminarayana  and  co-workers  at  Pennsylvania  State  University  at  the 
AIAA  2«th  Aerospace  Sciences  Meeting  in  Reno,  January  1986.  These 
procedures  are  currently  being  investigated  for  hydrodynamics  problems  at 
MARIN,  the  National  Maritime  Institute  of  the  Netherlands  by  H.  Raven  and  M. 
Hoekstra,  for  transonic  viscous  boattail  configurations  at  the  NASA  Langley 
Research  Center  with  R.  Wilmoth,  for  primitive  variable  formulations  with 
subsonic  viscous/ invisc id  interaction  at  ONR  with  T.C.  Tai  as  technical 
monitor,  at  the  Allison  Gas  Turbine  Division  of  General  Motors  by  D.  Reddy 
for  internal  flow  problems,  for  the  NASA  Lewis  Research  Center  under  J. 
Adamcyzk  and  B.  Anderson  for  internal  and  hypersonic  flows,  respectively. 
The  AFOSR  sponsored  work  has  also  been  referenced  in  numerous  papers  during 
the  past  year  and  very  similar  ideas  appear  in  the  work  of  P.  Bradshaw, 
Imperial  College  and  R.  Consteix  at  CERT,  France. 


It  has  been  approximately  five  years  since  the  formulations  considered 
here  were  first  proposed  by  the  present  investigators.  These  methods  have 
been  shown  to  be  accurate  and  efficient  procedures  for  two-dimensional 
steady  flows  when  combined  with  the  coupled  strongly  implicit,  conjugate 
gradient  and  global  relaxation  algorithms.  The  initial  applications  to  two- 
dimensional  unsteady  and  three-dimensional  steady  problems,  that  have  been 
reported  here,  have  further  established  the  utility  of  such  techniques. 
Until  such  time  as  a  fully  coupled,  time-dependent,  compressible  Navier- 
Stokes  solver  becomes  cost  and  computer  efficient  for  solving  general 
viscous  flow  problems,  the  procedures  discussed  herein  will  continue  to  be 
highly  competitive  for  a  significant  class  of  aerodynamic  configurations  and 


viscous/ inviscid  interactions 


